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1. 


INTRODUCTION 


Large  amounts  of  surface  gravity  data  now  exist  in 
organizations  such  as  the  Department  of  Defense  (DoD)  Gravity 
Library  located  at  the  Defense  Mapping  Agency  Aerospace  Center 
(DMAAC).  Present  holdings  are  expected  to  increase  substan¬ 
tially  in  the  future  as  data  from  other  sophisticated  sensor 
technologies  become  available.  Such  sensors  include  satellite 
radar  altimeters,  satellite-to-satellite  tracking  instrumenta¬ 
tion,  and  airborne  gravity  gradiometers .  For  example,  DMAAC 
holdings  of  GEOS-3  satellite  radar  altimeter  along-track  geoid 
heights  were  more  than  one  million  points  in  1978  even  after 
smoothing  to  a  1  sec  data  rate  (Ref.  3). 

Many  applications  require  an  accurate  estimation  of 
the  earth's  gravitational  field  from  the  existing  partial-earth 
coverage  gravity  information.  In  particular,  precise  values 
of  the  gravity  disturbance  vector  are  needed  for  the  compensa¬ 
tion  of  modern  inertial  guidance  and  navigation  systems.  An 
attractive  and  widely-known  technique  is  statistical  minimum- 
variance  estimation  (or  "least-squares  collocation").  This 
method  has  optimal  properties  with  respect  to  the  assumed 
models,  and  is  easily  adapted  to  a  wide  class  of  problems 
(Refs.  9  and  17).  The  major  inhibitor  to  the  widespread  appli¬ 
cation  of  optimal  gravity  data  processing  is  the  severe  com¬ 
puter  time  and  storage  demand  when  used  with  high-density 
gravity  surveys. 


The  concept  exploits  the  special  structure  of  gravity  model 
covariance  functions,  together  with  the  computational  efficiency 
of  the  Fast  Fourier  Transform  (Refs.  4  and  5),  to  achieve  a 

JL 

computational  complexity  of  order  Nlog2N."  The  result  is  an 

estimation  technique  which  is  suitable  for  processing  very 

large  gravity  databases.  In  contrast,  standard  methods  have  a 

3 

complexity  of  order  N  ,  and  the  best  competing  algorithm,  due 

2 

to  Levinson  (Refs.  1,  14,  16,  and  23),  is  of  order  N  complex¬ 
ity.  For  1000  data  points  such  methods  are  up  to  100  times 
more  expensive  than  the  present  technique  in  terms  of  computer 
time . 

The  basic  frequency  domain  approach  of  Generalized 
Wiener  Filtering  (Refs.  18  and  22)  assumes  planar  data  of  infi¬ 
nite  extent,  and  results  in  low  accuracy  when  truncated  to 
finite  data  sets.  An  extension  to  this  method  suitable  for 
finite  track  lengths  in  one-dimensional  estimation  was  de¬ 
scribed  in  Ref.  10.  The  error  properties  of  the  proposed 
method  were  analyzed  in  detail  and  shown  to  be  superior  to  the 
basic  technique,  while  retaining  the  same  order  of  complexity. 
This  method,  called  GEOFAST  for  "geodetic  fast  estimation", 
was  refined  and  implemented  using  simulated  along-track  grav¬ 
ity  data  with  good  results  (Ref.  21).  The  algorithm  has  now 
been  extended  to  two-dimensions  resulting  in  further  refinement 
and  modification,  while  retaining  its  accuracy  and  efficiency 
properties .  A  test  plan  was  developed  to  evaluate  the  perform¬ 
ance  of  GEOFAST  using  satellite  radar  altimeter  data  from 
GEOS-3,  and  has  been  completed  with  good  results. 

This  report  is  organized  into  four  chapters  and  two 
appendices.  Chapter  2  describes  the  two-dimensional  GEOFAST 


^Complexity  is  the  functional  dependence  of  the  number  of 
multiplications  (hence  computer  time)  on  the  number  of  data 
points,  N,  for  large  N. 
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algorithm  in  outline  form  with  reference  to  the  appendices  for 
details.  Appendix  A  is  a  brief  development  of  the  mathemati¬ 
cal  theory  underlying  the  one-dimensional  version  and  is  a 
condensation  of  material  in  Ref.  10.  Appendix  B  describes  the 
theory  necessary  to  extend  the  method  to  two  dimensions  includ¬ 
ing  a  definition  of  the  relevant  tensor  concepts.  The  results 
of  two  test  case  applications  using  GEOS-3  data  are  presented 
in  Chapter  3,  together  with  performance  evaluations.  The  last 
chapter,  Chapter  4,  contains  a  short  summary. 


2. 


THE  GEOFAST  ALGORITHM 


2.1  FORMULATION  OF  THE  PROBLEM 

The  GEOFAST  algorithm  was  developed  to  provide  an 
efficient  computational  solution  to  the  minimum-variance  esti¬ 
mation  equations 


In  geodetic  work,  this  formula  is  commonly  known  as  least-square 
collocation  (Ref.  9  and  17).  Here  z  is  a  data-vector  of  dimen¬ 
sion  N  and  C  its  auto-covariance  matrix  of  dimension  N*N. 
zz 

The  vector  of  estimates  x  is  of  dimension  M,  and  C  is  the 

cross-covariance  matrix  between  x  and  z  of  dimension  M*N.  The 

covariance  C  is  generally  modelled  by  an  analytic  function 

with  parameters  fitted  to  global  or  local  data,  and  includes 

measurement  noise.  The  cross-covariance  C  is  then  analyti- 

xz 

cally  derived  from  C  by  using  the  physical  relationship  be- 

Z  Z 

tween  the  estimated  quantities  x  and  the  observed  quantities  z 
(Refs.  11,  13,  and  24).  The  method  is  implicit  since  this  re¬ 
lationship  does  not  appear  explicitly  via  a  measurement  equation 
of  the  form  z  =  H  x.  Formally,  the  covariances  are  defined  by 

Czz=E!5zT}  ,  Cxz  =  Eli  zT)  (2.1-2) 

where  E  is  the  expectation  operator  and  both  vectors  are  assumed 
to  have  zero  means. 

Computer  solution  of  Eq .  2.1-1  by  standard  techniques 
involves  on  the  order  of  N~*  operations  (such  as  multiplications 


and  additions)  to  obtain  the  intermediate  vector 


S  =  C'zl  *  (2.1-3) 

and  on  the  order  of  MN  operations  to  yield  the  estimates 

x  =  Cxz  u  (2.1-4) 

This  limits  the  feasible  application  of  these  methods  to  data 
sets  of  a  few  hundred  points  at  most,  whereas  thousands  of 
measurements  are  often  available  from  such  sensors  as  satel¬ 
lite  altimeters. 

By  making  two  assumptions  and  using  frequency  domain 
techniques  this  workload  can  be  reduced  dramatically  to  the 
order  of  Nlog2N.  First,  let  the  data  points  z(6  ,<)>)  be  given 
on  a  rectangular  grid  in  the  6 , 0  plane  (Fig.  2.1-1).  Second, 
assume  that  the  covariance  between  two  data  points  is  a  func¬ 
tion  of  relative  position  only.  That  is,  the  covariance  func¬ 
tion  is  shift-invariant  (or  a  displacement  kernel) 

COV  [z(6^  ,<(>^) ,  z(02  -  ^  zz^®  2~® 1  ’  ^2-<^l^  (2*1-5) 

Note  that  this  assumption  includes,  but  is  not  restricted  to, 
the  case  of  isotropic  models.  Similar  assumptions  are  made 
for  the  cross-covariance  and  estimate  points  x(0,<t>),  namely, 

cov  ix(e1 ,4»1 ) ,  z(e2,<i>2)]  =  fxz(e2'ei’  ",2'<t,l)  (2.1-6) 

where  it  is  here  assumed  for  simplicity  that  the  estimates  are 
obtained  on  the  same  grid  as  the  data. 

The  foregoing  assumptions  are  not  unduly  restrictive 
in  practice  as  shown  by  the  application  to  GEOS-3  data  in  Chap¬ 
ter  3.  Data  given  on  a  limited  region  of  a  spherical  (or 


R -49247 


0 - ► 

COV  (P,Q)  =  COV  (P',Q')  =  fzz  (07~0y  02-0, ) 

Figure  2.1-1  Data  and  Covariance  Structure 

ellipsoidal)  surface  may  be  mapped  onto  a  plane  with  very  lit 
tie  distance  distortion  and  then  gridded.  The  shift  invari¬ 
ance  assumption  then  corresponds  to  the  common  model  of  sta¬ 
tionary  statistics.  The  estimate  grid  may  actually  be  any 
translation  (that  is  offset  in  altitude,  latitude,  and  lon¬ 
gitude)  of  the  data  grid.  With  a  fast  algorithm  it  is  often 
cheaper  to  compute  estimates  at  all  grid  locations  even  if 
only  a  subset  is  wanted. 

The  stated  assumptions  give  to  the  estimation  Eq . 
2.1-1  a  special  structure  which  is  exploited  by  the  GEOFAST 
algorithm.  The  natural  form  for  the  data  is  an  nj*n2  matrix 

Z  -  ,  z  j  k  —  z(0q  +  jA0,  $  Q+kA0 )  (2.1-7) 

where  j =0 , 1  ,  . . .  , n ^ - 1  and  k=0 , 1  ,  . . .  ,n2~l  corresponding  to  the 
data  grid  (Fig.  2.1-1)  with  origin  (0q,0q)  and  mesh  size  (A0  , 
A0 ) .  The  total  number  of  grid  points  is  N  =  n-^-  The  data 
matrix  Z  is  put  in  vector  form  z  by  listing  its  elements  in 
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row-by-row  order.  With  this  convention  the  N*N  covariance 
matrix  Czz  assumes  a  block  structure  (see  Appendix  B) 


in  which  there  are  n^xn^  blocks  and  each  block  is  of  dimension 
n2xn2*  The  block  matrix  is  the  cross-covariance  between  the 
elements  in  two  rows  of  Z  which  are  k  rows  apart.  This  depends 
only  on  k  by  shif t- invariance .  By  the  same  argument  each  block 
matrix  is  of  the  form 


The  matrix  C  is  said  to  be  block  Toeplitz  with 
Toeplitz  blocks,  or  simply  block  Toeplitz.  Note  that  Czz  is 
symmetric  and  is  completely  determined  by  its  first  block  row 
(or  block  column),  whereas  the  blocks  (except  Tq)  are  not  sym¬ 
metric  and  are  determined  by  their  first  row  and  column.  Alto¬ 
gether,  the  definition  of  Czz  requires  less  than  2n^n2  =  2N  pa¬ 
rameters.  The  structure  of  C  is  the  same  except  that  symmetry 

X  z 

does  not  hold  in  general,  so  that  approximately  4N  parameters 
are  needed  to  define  it.  In  summary,  the  geodetic  estimation 
problem  is  equivalent  to  the  inversion  and  multiplication  of 
block  Toeplitz  matrices. 


2.2  OUTLINE  OF  THE  ALGORITHM 


The  GEOFAST  algorithm  achieves  its  efficiency  by  trans¬ 
forming  the  estimation  equations  into  the  frequency  domain  where 
an  accurate  approximation  can  be  made  to  reduce  the  workload. 

The  transformation  is  accomplished  with  the  Fast  Fourier  Trans¬ 
form  (FFT)  which  requires  on  the  order  of  Nlog2N  operations. 

The  transformed  covariance  matrices  are  closely  approximated 
by  (block)  banded  matrices  resulting  in  a  solution  workload  of 
order  N.  A  simple  trade-off  between  approximation  accuracy  and 
computer  workload  is  controllable  by  choice  of  the  bandwidth. 

A  simplified  flow  chart  of  the  algorithm  is  shown  in 

Fig.  2.2-1.  Input  to  the  computation  consists  of  the  n^xn2  data 

matrix  Z,  together  with  the  definition  of  the  covariance  matrices 

C  ,  C  .  As  seen  in  Eqs.  2.1-8  and  2.1-9,  each  covariance  can 
xz  zz  n 

be  specified  by  4n^2  parameters  which  are  naturally  arranged 
as  a  2n^x2n2  correlation  matrix.  Specifically,  let  the  Toeplitz 
matrix  T  in  Eq.  2.1-9  be  represented  by  the  column  vector 

°’  t-(n2-l)’  •••’  t-l*  t0>  tl*  •••’  tn2-l]T 

(2.2-1) 


of  dimension  2n2 .  Then  each  block  matrix  T^  in  Eq.  2.1-8  can 
be  represented  by  a  vector  t_^ ,  and  the  matrix  Czz  by  the  2n^x2n2 
matrix 


zz 
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(2.2-2) 


■’'Superscript  T  is  used  throughout  this  report  to  denote  the 
transpose  operator.  When  not  a  superscript,  T  is  used  to  de¬ 
note  a  Toeplitz  matrix  without  causing  ambiguity. 
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Figure  2.2-1  Two-Dimensional  GEOFAST  Algorithm 


For  the  symmetric  case  (Czz>,  only  half  of  the  matrix  T^z  is 
required  with  t_k  h  t^.  Total  computer  storage  for  input  is 
of  order  N  =  n-j^- 

2.2.1  Transformation  to  the  Frequency  Domain 

The  first  stage  of  the  algorithm  is  the  transformation 
of  the  inputs  to  the  frequency  domain.  The  FFT  is  simply  an 
efficient  implementation  of  the  complex  Discrete  Fourier  Trans- 
form  (DFT)  which  is  defined  in  two  dimensions  by 


1 

Vni  n2 


n2 

X  zjk  exe 

k=0 


I 


-  2n  i 


(2.2-3) 


Since  the  estimation  equations,  Eq .  2.1-1,  are  in  terms  of  the 
data  vector  z  it  is  convenient  to  express  Eq .  2.2-3  in  the 
equivalent  matrix  form, 

*The  symbol  i  denotes  the  imaginary  unit  fT. 
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(2.2-4) 


z'  =  F  z 

The  N*N  complex  two-dimensional  DFT  matrix  F  has  a  special 
structure  which  is  discussed  further  in  Appendix  B.  The  vec¬ 
tor  z'  in  Eq .  2.2-4  is  the  complex  representation  of  the  data 
vector  z  in  a  finite  Fourier  series  of  complex  exponentials. 

Since  real  quantities  are  more  convenient  to  deal  with  in  a 
computer,  a  further  transformation  to  a  sine  and  cosine  series 
is  performed  by  combining  real  and  imaginary  parts.  This  cor¬ 
responds  to  the  matrix  operations  (redefining  z' ) 

z'  =  H  F  z  (2.2-5) 

where  H  is  a  very  sparse  complex  N*N  matrix.  The  product  HF 
is,  of  course,  a  real  matrix. 

Finally,  the  accuracy  of  the  banded  approximation, 
which  is  essential  for  the  success  of  the  algorithm,  is  criti¬ 
cally  enhanced  by  the  introduction  of  a  data  window.  In  matrix 
form  this  becomes  (redefining  z'  again) 

z'  =  H  F  W  z  (2.2-6) 

where  W  is  a  real  N*N  diagonal  matrix.  That  is,  each  data  value 
z^  is  multiplied  by  a  scalar  weight  w^  before  applying  the  DFT. 
This  windowing  compensates  for  the  finite  extent  of  the  data  grid 
in  a  manner  explained  in  Appendix  A.  A  representative  window 
function  in  two  dimensions  would  be  a  gaussian  probability  den¬ 
sity  which  is  peaked  at  the  center  and  tapered  toward  the  edges 
of  the  data  grid.  In  GEOFAST,  an  optimal  one-dimensional  window 
w°  derived  by  Kaiser  (Ref.  15)  is  used  to  form  a  product  window 
Wj  k  =  (Wj)(w£)  (Fig.  2.2-2).  The  combined  data  transformation 
in  Eq .  2.2-6  can  be  implemented  in  order  Nlog2N  operations 
since  multiplication  by  H  and  W  are  both  order  N. 
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Figure  2.2-2  Product  Window  Definition 

Corresponding  to  the  frequency  domain  data  vector  z’ 
(consisting  of  Fourier  coefficients)  is  the  auto-covariance 
matrix 


CzZ  =  E{z'(z’)T}  =  (H  F  W)  Czz  (H  F  W)T  (2.2-7) 

Similarly,  introducing  the  transformed  estimates 

x'  = (h  F  w)x  (2.2-8) 

leads  to  the  cross-covariance  matrix 

Cxz  =  eU'(z')T}  =  (H  F  W)  Cxz  (H  F  W)T  (2.2-9) 

A  straightforward  row  and  column  implementation  of  Eqs.  2.2-7 

2 

and  2.2-9  would  result  in  an  order  N  log2N  algorithm.  This  is 
avoided  by  replacing  C and  C'  by  their  banded  approximations 
and  calculating  only  those  elements  within  the  retained  bands. 
If  the  number  of  bands  in  the  one  dimensional  window  is  mD  this 


2 

reduces  the  workload  to  the  order  of  mgNlog2N,  where  for  typi¬ 
cal  applications  mg  is  small  and  independent  of  N.  The  scheme 
used  to  compute  C'  and  C'  bandwise  is  derived  in  Appendices 

XZ  Z  Z 

A  and  B,  and  is  an  essential  part  of  the  GEOFAST  algorithm. 
2.2.2  Solution  of  the  Estimation  Equations 

The  preceding  transformations  constitute  the  first 

stage  of  the  calculation  (Fig.  2.2-1).  Note  that  storage  re- 

2 

quirements  for  the  covariances  are  of  order  mgN .  The  second 
stage  consists  of  the  solution  of  the  estimation  equations  in 
the  frequency  domain.  The  estimation  equations  in  the  trans¬ 
formed  variables  become  simply 

*’  *  C^Czz)'1  5'  (2.2-10) 

coupled  with  the  inverse  of  Eq .  2.2-8 

x  =  (H  F  W)'1  x'  (2.2-11) 

It  is  easily  verified  that  the  solution  of  Eqs.  2.2-10  and 
2.2-11  is  identical  to  the  solution  of  Eq .  2.2-1.  As  before, 
the  equations  are  solved  in  two  steps 

u'  =  (C-z)-:  z'  (2.2-12) 

x'  =  C'  u'  (2.2-13) 

—  XZ  — 

Since  C'  is  banded,  the  multiplication  in  Eq .  2.2-13  requires 
x  z  2 

only  on  the  order  of  m^N  operations. 


However,  analysis  of  the  solution  of  Eq .  2.2-12  is 
more  complicated.  It  is  shown  in  Appendix  B  that  the  block 


Toeplitz  structure  of  C  in  Eq.  2.1-8  gives  rise  to  a  corre- 

2  2 

sponding  block-banded  structure  for  C22  in  which  the  blocks 
are  also  banded  (Fig.  2.2-3).  There  is  an  essential  differ¬ 
ence  between  this  structure  and  the  simple  band  structure  of  a 
diagonal  block  which  is  the  one-dimensional  analogue  of  C^2 . 

In  the  one-dimensional  problem, the  symmetric  band  matrix  C' 

2  22 
can  be  inverted  in  order  mgn  operations  by  the  standard  Cholesky 

technique  (Refs.  6  and  12),  where  n  is  the  matrix  dimension  and 
mg  the  number  of  superdiagonal  bands.  When  the  Cholesky  tech¬ 
nique  is  applied  to  the  block-banded  structure  of  the  two- 
dimensional  problem,  the  internal  banding  within  the  blocks  is 
lost  during  the  calculation.  This  difficulty  is  known  as 
"fill-in".  The  result  is  an  effective  bandwidth  of  ,  where 

mg  is  now  the  number  of  blockbands  of  size  n2*n2»  and  a  work¬ 
load  of  order  (mgn2  ^(n-j^ )  •  Assuming  n^  =  n2  gives  an  order 
of  effort  proportional  to  n^  =  rather  than  N. 


R-49251 


Figure  2.2-3  Block  Banded  Structure 

An  order  N  algorithm  for  the  solution  of  Eq .  2.2-12 
can  be  developed  by  using  an  iterative  method  based  on  a  special 


class  of  covariance  matrices.  This  special  type  of  covariance 
structure  is  termed  separable  and  is  characterized  by  a  covar¬ 
iance  function  in  Eq.  2.1-5  of  the  form 

f ( 0 2 =  ^1^®2”^1^  ^2^2~^1  ^  (2.2-14) 

This  class  of  functions  is  not  general  enough  to  include  the 
statistical  models  required  for  gravimetric  applications,  but 
it  can  serve  as  an  approximation  for  use  in  an  iterative  method. 
It  is  shown  in  Appendix  B  that  for  a  separable  covariance  ma¬ 
trix  Eq.  2.2-12  "factors"  into  two  one-dimensional  matrix  equa¬ 
tions  of  the  form 


V'  =  (C^)'1(ZI )T  (2.2-15) 

U'  =  (C])'1(V')T  (2.2-16) 

where  Cj ,  correspond  to  f^ ,  f2  and  are  square  band  matrices 
of  dimension  n^  ,  1^2-  The  matrices  Z'  and  U'  are  both  of  dimen¬ 
sion  n^*n2>  and  are  equivalent  to  the  vectors  z'  and  u'  when 
their  elements  are  listed  in  row-by-row  order.  If  the  band¬ 
width  of  the  matrix  Cl,  is  mg,  then  Eq .  2.2-15  can  be  solved  by 

the  Cholesky  method  in  order  mgn2(mg+n^)  operations.  Assuming 

nl  =  n2  ’  t*ie  total  workload  for  both  equations  is  thus  of  order 
nigN  in  the  separable  case. 

The  iterative  solution  of  Eq .  2.2-12  in  the  non- 
separable  case  proceeds  as  follows.  Choose  the  matrix  D'  as  a 
separable  approximation  to  C^z  (which  may  be  done  in  a  natural 
way)  and  define  the  matrix  E'  by 

C'  =  D'  +  E'  (2.2-17) 

zz 
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The  matrix  E'  will  be  block-banded  since  both  C'  and  D'  are. 

zz 

The  estimation  equation  can  now  be  written  implicitly 


u'  =  (D ’ )_1  [ z 1  -  E ' u ' ] 


(2.2-18) 


or  in  recursive  form 

u ' (k+1 )  =  (D')'1  [z'  -  E'u'(k)]  (2.2-19) 

where  k  indicates  the  iteration  number  and  u'(0)  is  any  con¬ 
venient  initial  solution. 

Since  both  the  inversion  and  multiplication  in  Eq. 
2.2-19  require  on  the  order  of  N  operations,  so  does  one  step 
of  the  iteration.  If  the  number  of  steps,  m  ,  necessary  for 

b 

convergence  of  the  iteration  is  independent  of  N  and  small 
enough,  the  total  workload  to  solve  Eq .  2.2-12  is  of  order  N. 

It  can  be  shown  that  the  rate  of  convergence  depends  on  the 
norm  of  the  matrix  ( D ’ )  ^E'  which  will  be  small  if  the  separa¬ 
ble  approximation  is  a  good  one.  Furthermore,  a  simple  modi¬ 
fication  to  the  iteration  guarantees  convergence  for  any  posi¬ 
tive  definite  matrix  D'.  The  convergence  properties  of  the 
iterative  technique,  Eq.  2.2-19,  and  the  appropriate  choice 
of  the  approximating  matrix  D'  are  discussed  in  Appendix  B. 

The  solution  of  the  estimation  equations  in  the  fre¬ 
quency  domain  completes  the  second  stage  of  the  GEOFAST  algo¬ 
rithm  (Fig.  2.2-1).  The  third  stage  consists  of  the  inverse 
transformation  back  to  the  space  domain  (Eq.  2.2-11).  Since  W 
is  a  diagonal  matrix,  its  inverse  is  an  order  N  operation.  The 
complex  DFT  matrix  F  can  be  inverted  in  order  Nlog2N  operations 
by  means  of  the  inverse  FFT ,  while  the  inverse  of  the  matrix  H 
can  be  obtained  analytically  and  has  the  same  sparse  structure 
as  H  itself.  Thus  the  final  stage  has  a  workload  dominated  by 
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The  computer  time  and  storage  requirements  for  the 
GEOFAST  algorithm  are  summarized  in  Table  2.2-1.  The  numbers 
given  are  asymptotic  estimates,  valid  for  large  N,  and  do  not 
include  multiplicative  factors  independent  of  N  except  for  mB 
and  mg .  For  repeated  application  of  the  algorithm  with  the 
same  covariance  model,  the  matrices  C^z  and  C^z  can  be  stored, 
reducing  the  time  requirement  to  the  order  of  Nlog2N  for  each 
application  to  a  new  data  set. 


TABLE  2.2-1 

COMPUTER  TIME  AND  STORAGE  REQUIREMENTS 


ALGORITHM 

STAGE 

ASYMPTOTIC  ORDER  OF  MAGNITUDE 

TIME 

STORAGE 

1: 

Input 

Transforms 

2 

mgNlog2N 

■4n 

2: 

Estimation 

Equations 

msmBN 

3: 

Output 

Transform 

Nlog2N 

N 

Full  Solution 

m^Nlog2N 

mBN 

N 

=  Total  number  of  data  points 

nig  =  Bandwidth  in  one  dimension 
mg  =  Number  of  steps  for  convergence 
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APPLICATION  TO  GEOS -3  DATA 


In  order  to  evaluate  the  GEOFAST  technique  a  Test 
Plan  was  developed  in  coordination  with  DMAAC.  The  objectives 
of  the  plan  were  to  test  the  effectiveness  of  GEOFAST  when 
applied  to  gravity  survey  data  and  to  compare  accuracy  and 
efficiency  with  existing  DMAAC  methods.  It  was  decided  that 
GEOS-3  satellite  altimeter  data  would  be  a  natural  choice  for 
the  test  plan  since  a  large  data  set  existed,  and  computer 
processing  time  is  a  critical  parameter  in  the  ability  to 
utilize  such  information  in  an  optimal  manner. 

A  20  deg  by  20  deg  area  in  the  Atlantic  was  chosen 
from  which  to  process  GEOS-3  data.  The  test  area,  bounded  by 
17  deg  and  37  deg  N  latitude  and  285  deg  and  305  deg  E  longi¬ 
tude  (Fig.  3-1)  includes  Bermuda  and  Puerto  Rico.  Much  of  the 
GEOS-3  calibration  area  is  also  included  so  that  the  track 
density  is  relatively  high.  The  geoid  heights  derived  from 
the  altimeter  measurements  were  then  used  to  estimate  mean 
gravity  anomalies.  The  implementation  of  the  test  plan  re¬ 
sulted  in  the  two  test  cases  described  below. 

The  computer  program  used  to  generate  the  test  case 
results  differed  in  some  aspects  from  the  algorithm  described 
in  Chapter  2.  The  test  program  utilized  a  preliminary  version 
of  GEOFAST  in  which  the  iterative  scheme  of  Section  2.2.2  was 
not  operational.  An  alternative  approach  was  implemented  re¬ 
quiring  increased  computer  storage  to  achieve  comparable  accu¬ 
racy  and  efficiency.  A  full  implementation  of  the  GEOFAST 
algorithm  is  currently  being  developed,  and  will  be  subjected 
to  the  same  test  plan  when  available. 
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Figure  3-1 


GEQS-3  Test  Area 


3.1  TEST  CASE  NO.  1 

Basic  data  for  this  test  case  were  the  1  deg  gridded 
geoid  heights  derived  by  DMAAC  (Ref.  3)  with  an  estimated  accu¬ 
racy  of  approximately  1  m.  Mean  gravity  anomaly  estimates 
were  obtained  on  the  same  1  deg  grid  for  the  central  10  deg  by 
10  deg  square  of  the  test  area.  Three  different  algorithms 
were  used  to  obtain  the  estimates:  GEOFAST,  Spatial  Colloca¬ 
tion,  and  Inverse  Stokes  (Ref.  9),  the  latter  two  being  tech¬ 
niques  in  use  at  DMAAC.  All  algorithms  used  exactly  the  same 
input  data  so  that  comparisons  could  be  made,  and  in  the  case 
of  GEOFAST  and  Spatial  Collocation  the  same  DMAAC-deve loped 
covariances  were  used. 

Two  transformations  were  applied  to  the  basic  data  to 
produce  the  input  data  for  the  algorithms.  First,  a  global 
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spherical  harmonic  expansion  of  the  geoidal  undulation  to  de¬ 
gree  and  order  12  was  subtracted  from  the  GEOS-3  geoid  heights. 
These  are  termed  residual  GEOS-3  geoid  heights  and  are  contour- 
plotted  in  Fig.  3.1-l(a).  Second,  a  "template"  approach  was 
used  to  generate  the  estimates.  That  is,  the  10  deg  by  10  deg 
square  for  which  the  mean  gravity  anomalies  are  desired  was 
subdivided  into  four  5  deg  by  5  deg  subregions  (Fig.  3.1-2). 

Within  each  subregion,  only  the  surrounding  15  deg  by  15  deg 

J-. 

square  of  data  was  used  in  forming  the  estimates  (Fig.  3.1-3). 

For  each  15  deg  square,  the  mean  residual  geoid  height  was  com¬ 
puted  and  subtracted  from  the  data  to  remove  bias.  The  resid¬ 
ual  mean  gravity  anomaly  estimates  were  generated  independently 
by  subregion  and  combined  to  yield  estimates  over  the  10  deg 
region.  To  these  residual  estimates,  (Fig.  3.1-l(b)),  the 
spherical  harmonic  expansion  of  the  (point)  gravity  anomaly  to 
degree  and  order  12  was  added  to  obtain  the  final  mean  gravity 
anomaly  estimates. 

In  the  statistical  methods  (GEOFAST  and  Spatial  Col¬ 
location)  ,  the  covariances  for  the  residual  geoid  heights  and 
residual  mean  gravity  anomalies  were  obtained  from  the  Tscherning- 
Rapp  model  No.  4  (Ref.  24)  with  terms  to  degree  and  order  12 
removed.  The  mean  gravity  anomalies  were  modeled  as  point 
anomalies  at  10.4  km  altitude.  (See  Ref.  24  for  a  justifica¬ 
tion  of  this  approximation . )  The  resulting  covariances  are 
displayed  as  a  function  of  spherical  distance  in  Fig.  3.1-4. 

For  an  independent  reference,  the  estimated  mean  grav¬ 
ity  anomalies  were  compared  with  values  from  the  DoD  Gravity 
Library  obtained  by  surface  gravity  surveys.  This  data  base 
is  on  a  5  min  grid  and  was  averaged  over  ]  deg  squares  to  yield 
comparison  values.  Statistics  of  the  comparison  are  displayed 

*For  the  Inverse  Stokes  method,  only  a  9  deg  by  9  deg  square 
was  used. 
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LONGITUDE  (deg  E) 


(b)  GEOFAST  Residual  Mean  Gravity  Anomaly  Estimates  (mgal) 


Figure  3.1-1 


Residual  GEOS-3  Geoid  Heights  and  GEOFAST 
Residual  Mean  Gravity  Anomaly  Estimates 
for  Test  Case  No.  1 


Figure  3.1-2  Subdivision  of  the  Estimation  Region 
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COLLOCATION 


NOTE:  IN  ALL  CASES  15°  x  15°  RESIDUAL  GEOID  MEAN  SUBTRACTED  FROM  DATA 


Figure  3.1-3  Estimation  Templates 
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(a)  Residual  Geoid  Height  Auto-Covariance 
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(b)  Residua]  Mean  Gravity  Anomaly/Geoid 
Height  Cross-Covariance 


Figure  3 . 1  - A 


Model  Covariances 


in  Table  3.1-1.  The  mean  discrepancy  is  less  than  1  mgal  and 
the  standard  deviation  about  4  mgal.  The  distribution  of  the 
differences  closely  approximates  a  normal  distribution.  It 
should  be  emphasized  that  the  observed  differences  include 
four  separate  effects: 

•  Errors  in  GEOS-3  geoid  height  measure¬ 
ments 

•  Incorrect  statistical  modeling 

•  Errors  due  to  finite  data  extent  and 
discretization 

•  Errors  in  surface  gravity  survey  data 


Any  errors  introduced  by  approximations  in  the  GEOFAST  algo¬ 
rithm  itself  are  too  small  to  be  seen  in  this  comparison  (See 
Ref.  21). 


TABLE  3.1-1 

STATISTICS  OF  THE  DIFFERENCE  BETWEEN  GEOFAST  MEAN 
GRAVITY  ANOMALY  ESTIMATES  AND  SURFACE  REFERENCE  DATA 


REGION 
(Fig.  3.1-2) 

MEAN 
(mgal ) 

STANDARD 
DEVIATION 
(mgal ) 

MAXIMUM 

MAGNITUDE 

(mgal) 

10°  x  10 

0 

0.6 

4.1 

11.7 

A 

1.3 

2.8 

6.7 

5°  x  5° 

B 

-3.1 

3.1 

-8.5 

C 

3.8 

3.9 

8.2 

D 

0.4 

3.1 

11.7 

Comparisons  were  also  made  at  DMAAC  with  estimates 
from  the  other  two  methods  in  use.  Table  3.1-2  summarizes 
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TABLE  3.1-2 

COMPARISONS  BETWEEN  ALGORITHMS 


RMS  DIFFERENCE 
(mgal ) 

GEOFAST 

SPACF  ! 

COLLOCATION 

INVERSE 

STOKES 

SURFACE 

DATA 

GEOFAST 

- 

00 

3.6 

mam 

SPACE 

COLLOCATION 

- 

3.1 

3.5 

INVERSE 

STOKES 

- 

3.5 

these  results,  including  a  comparison  between  each  algorithm 
and  surface  data.  Space  Collocation  and  Inverse  Stokes  both 
agreed  slightly  better  with  the  surface  reference  than  GEOFAST 
(3.5  versus  4.2  mgal ) ,  but  differed  from  each  other  by  almost 
the  same  amount  (3.1  mgal).  GEOFAST  differed  from  Space  Col¬ 
location  by  less  than  2  mgal.  Although  these  two  methods  are 
based  on  the  same  assumptions,  there  was  one  known  discrepancy 
in  this  test  case.  The  GEOFAST  covariance  construction  ignored 
the  fact  that  the  E-W  distance  between  grid  points  decreases 
with  increasing  latitude,  instead  using  an  average  distance 
valid  at  27  deg  N  latitude.  This  introduces  an  error  of  5-10% 
which  could  account  for  up  to  1  mgal  rms  difference.  This 
source  of  error  was  avoided  on  the  second  test  case  by  using  a 
different  construction  (Section  3.2). 

Overall,  the  results  verified  the  accuracy  of  the 
GEOFAST  estimates  and  the  suitability  of  the  method  for  possi¬ 
ble  use  with  GEOS-3  geoid  height  data.  The  relatively  small 
size  of  this  first  test  case  (N=225)  did  not  admit  meaningful 
timing  and  efficiency  comparisons  to  be  made. 
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3.2  TEST  CASE  NO.  2 


The  objective  of  this  test  case  was  to  demonstrate 
the  ability  of  GEOFAST  to  process  high  density  along-track 
data.  Accordingly,  the  database  for  the  second  test  case  was 
the  GEOS-3  adjusted  along-track  point  geoid  height  library 
(Ref.  3).  This  data  has  been  smoothed  along-track  to  a  1  sec 
(7  km)  spacing,  and  adjusted  to  minimize  track  crossing  dis¬ 
agreement.  Data  contained  in  the  same  20  deg  by  20  deg  square 
used  in  Test  Case  No.  1  were  obtained  from  DMAAC  and  amounted 
to  about  90,000  points.  This  was  the  base  data  from  which  the 
1  deg  grid  values  were  obtained  for  the  previous  case.  Resid¬ 
ual  geoid  heights  were  obtained  as  before  by  subtracting  a 
spherical  harmonic  expansion  to  degree  and  order  12.  This 
reference  geoid  was  computed  on  a  30  min  grid  and  all  data 
points  within  a  single  grid  square  were  compensated  identi¬ 
cally  . 


A  transformation  of  coordinates  was  made  in  order  to 
obtain  gridded  data  for  the  GEOFAST  algorithm  in  which  the 
distance  between  points  is  constant  in  both  directions.  A 
Lambert  projection  (Ref.  7)  was  used  to  map  a  spherical  cap  of 
15  deg  radius  centered  on  the  test  area,  onto  a  circular  disc 
in  the  plane.  The  spherical  coordinates  of  latitude  (0)  and 
longitude  (A)  were  thereby  mapped  into  cartesian  coordinates 
(x,y)  in  such  a  way  that  distance  distortion  was  minimized. 

The  equations  of  this  transformation  are 


X 

sin 

% 

0 

-cos 

*0 

cos 

*0 

sin  Kq 

0 

cos  0 

COS 

A 

y 
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0 

1 
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-sin 
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cos  Aq 

0 

COS  0 

sin 
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z 

cos 

*0 

0 
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0 

0 

1 

sin 

0 

(3.2-1 ) 
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x  =  py 


(3.2-2) 


L'M'jilJI  UWIP 


.  y  =  -px 


where  <J>q  =  27  deg  and  Aq  =  295  deg,  and  result  in  less  than 
0.5%  distance  distortion  for  a  15  deg  cap.  The  mapping  is 
displayed  graphically  in  Fig.  3.2-1,  where  x  and  y  are  meas¬ 
ured  east  and  north,  in  degrees  of  great  circle  distance,  from 
the  origin  at  27  deg  N  latitude  and  295  deg  E  longitude.  The 
data  were  then  further  restricted  to  a  15  deg  by  15  deg  square 
in  the  plane  whose  exact  location  was  determined  by  the  avail¬ 
ability  of  surface  reference  data.  Figure  3.2-2  shows  the 
reduced  region  with  the  origin  shifted  to  coincide  with  x=-8, 
y=-6  in  Fig.  3.2-1.  The  mean  residual  geoid  height  over  this 
region  was  subtracted  from  the  data  to  remove  bias. 

The  second  test  case  was  split  into  two  subcases  for 
purposes  of  timing  and  accuracy  comparisons.  These  were 

•  Case  2A:  1  deg  by  1  deg  grid  (N=225) 

•  Case  2B:  30  min  by  30  min  grid  (N=900) 


where  N  represents  the  total  number  of  grid  points  in  the  15 
deg  square.  The  gridded  geoid  heights  were  obtained  by  simple 
averaging  with  respect  to  x,y  coordinates. 

The  gridded  data  thus  obtained  was  processed  by  GEOFAST 
(in  a  single  computation)  to  produce  residual  mean  gravity 
anomalies  at  the  same  grid  points  as  the  residual  geoid  heights. 
The  model  covariance  functions  were  the  same  as  those  used  in 
the  first  test  case  (Fig.  3.1-4).  For  surface  reference  data 
the  5  min  gravity  library  anomalies  were  reduced  by  the  spher¬ 
ical  harmonic  model  to  degree  and  order  twelve,  and  then  mapped 
onto  a  plane  and  gridded  in  exactly  the  same  manner  as  the 
GEOS-3  data.  Although  in  practice  the  mean  gravity  anomaly 
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Figure  3.2-2  Reduced  Estimation  Region 
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estiratates  would  be  mapped  back  to  latitude  and  longitude  co¬ 
ordinates,  this  was  not  required  for  the  present  purposes. 

The  residual  GEOS-3  geoid  heights  for  both  grid  spac- 
ings  (1  deg  and  30  min)  are  shown  in  Fig.  3.2-3,  and  the  GEOFAST 
residual  mean  gravity  anomaly  estimates  are  shown  in  Fig.  3.2-4. 
The  increased  resolution  at  the  higher  density  is  apparent, 
particularly  in  the  region  near  Bermuda  and  in  the  Bahamas. 

The  estimates  were  differenced  with  the  surface  reference  re¬ 
sidual  anomalies  and  the  statistics  of  this  comparison  are 
given  in  Tables  3.2-1  and  3.2-2,  Three  concentric  square  re¬ 
gions  were  considered:  the  central  5  deg  square,  the  central 
10  deg  square,  and  the  15  deg  square.  As  expected  the  accuracy 
of  the  estimates  decreased  toward  the  edge  of  the  finite  data 
base.  The  mean  difference  in  all  cases  was  less  than  1  mgal . 

The  difference  standard  deviation  of  3.2  mgal  for  the 
5  deg  region  with  the  1  deg  grid  agreed  closely  with  the  result 
from  Test  Case  No.  1.  This  increased  to  only  4.1  mgal  with  the 
30  min  grid,  whereas  theory  predicts  a  >/N"  dependence  for  a  fac¬ 
tor  of  two  increase.  Over  the  15  deg  region  the  values  were 
7.5  mgal  and  10.7  mgal,  respectively,  which  is  an  increase  of 
about  40%.  Maximum  differences  were  dominated  by  the  presence 
of  Bermuda  which  is  not  adequately  modeled  in  the  covariance 
structure.  Even  so,  this  singularity  is  estimated  with  much 
sharper  definition  at  the  higher  density  grid.  Also  noteworthy 
is  the  good  accuracy  obtained  with  gridding  by  simple  averaging, 
although  this  property  depends  on  the  smoothness  of  the  data. 

The  geographical  distribution  of  the  differences  between  sur¬ 
face  data  and  GEOFAST  estimates  was  also  analyzed.  Some  of 
the  difference,  although  clearly  not  all,  at  the  lower  and 
upper  boundaries  in  Case  2B  appears  to  be  due  to  GEOFAST  edge 
error  ( See  Ref .  21 ) . 


TABLE  3.2-1 

DIFFERENCES  BETWEEN  GEOFAST  MEAN  GRAVITY  ANOMALY  ESTIMATES 
AND  SURFACE  REFERENCE  DATA  FOR  CASE  2A  (1  deg  GRID) 


REGION 

(CENTERED) 

MEAN 
(mgal ) 

STANDARD 
DEVIATION 
(mgal ) 

MAXIMUM 
MAGNITUDE 
(mgal ) 

5°x5° 

-0.3 

3.2 

6.3 

10°xl0° 

0.2 

5.0 

28.6* 

15°x 15° 

0.8 

7.5 

28.6* 

Maximum  occurs  at  the  Bermuda  rise. 


TABLE  3.2-2 

DIFFERENCES  BETWEEN  GEOFAST  MEAN  GRAVITY  ANOMALY  ESTIMATES 
AND  SURFACE  REFERENCE  DATA  FOR  CASE  2B  (30  mTn  GRID) 


REGION 

(CENTERED) 

MEAN 
(mgal ) 

STANDARD 
DEVIATION 
(mgal ) 

MAXIMUM 
MAGNITUDE 
(mgal ) 

5°x5° 

-0.4 

4.1 

15.6 

10°xl0° 

0.1 

8.8 

JL 

136. 

15°xl5° 

0.7 

10.7 

136.* 

Maximum  occurs  at  the  Bermuda  rise. 


With  two  different  grids,  computer  time  comparisons  can 
be  made  (Table  3.2-3).  As  indicated,  the  computer  time  used  in 
pre-processing  to  obtain  gridded  data  was  modest  and  nearly  the 
same  for  both  grids.  The  total  Lime  used  by  the  GEOFAST  algo¬ 
rithm  increased  by  a  factor  of  about  six  when  the  number  of 
grid  points  increased  by  a  factor  of  four.  The  analysis  in 
Section  2  predicts  an  asymptotic  Nlog?N  dependence  which  would 
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TABLE  3.2-3 

GEOFAST  ALGORITHM  COMPUTATIONAL  EFFICIENCY 


STAGE 

TIME  ON  IBM  370  (min) 

CASE  2A 
(N=225 ) 

CASE  2B 
(N=900 ) 

PRE-PROCESSING 

<2 

<2 

TRANSFORM  COVARIANCES 

~1 

~5 

COMPUTE  ESTIMATES 

~1 

~7 

TOTAL  GEOFAST 

~2 

~12 

result  in  a  ratio  of  about  five.  Considering  the  relatively 
small  size  of  N  this  is  acceptable  agreement.  The  required 
time  to  compute  900  estimates  was  12  min  on  an  IBM  370/3033 
computer,  and  would  correspond  to  less  than  2  min  on  a  CDC 
7600.  The  results  of  this  test  justify  the  conclusion  that 
the  GEOFAST  algorithm  has  the  ability  to  process  high  density 
data  with  modest  pre-processing,  and  to  produce  high  resolu¬ 
tion  estimates  of  good  quality  with  considerable  efficiency. 


-'Estimated  time  on  a  CDC  7600  computer  is  based  on  benchmark 
testing  conducted  by  the  author. 
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4. 


SUMMARY 


An  algorithm  has  been  presented  which  allows  the 
rapid  calculation  of  minimum-variance  estimates  from  high- 
density  gravity  data.  The  GEOFAST  method  uses  frequency  do¬ 
main  techniques  to  achieve  a  complexity  of  order  Nlog2N  where 
N  is  the  number  of  data  points.  This  efficiency  is  obtained 
by  exploiting  two  structures  commonly  found  in  gravimetric 
applications : 

•  Input  data  that  is  regularly  spaced  on  a 
cartesian  grid  (after  pre-processing) 

•  Statistical  covariance  models  that  are 
stationary  (or  shift-invariant). 

The  resulting  algorithm  has  been  implemented  in  two-dimensions 
and  is  well  suited  to  the  processing  of  very  large  gravity 
data  bases.  A  balance  between  accuracy  and  computer  time  may 
be  controlled  by  the  choice  of  design  parameters. 

The  GEOFAST  algorithm  has  been  used  to  estimate  mean 
gravity  anomalies  from  about  90,000  GEOS-3  geoid  height  values 
suitably  gridded  to  900  points  on  a  15  deg  square  at  a  30  min 
resolution  (Fig.  4-1).  Agreement  with  surface  gravity  data 
was  relatively  good  and  computer  running  time  verified  the 
theoretical  bounds.  It  is  possible  to  extend  the  algorithm 
further  so  that  the  two  basic  assumptions  need  hold  only  approx 
imately,  thus  considerably  enhancing  the  range  of  application. 


43 


«-*77?4 


igure  4-1 


GEOFAST  Residual  ^Mean  Gravity  Anomaly 
Estimates  at  30  min  Resolution  (mgal) 


APPENDIX  A 

GEOFAST:  ONE-DIMENSIONAL  THEORY 


This  appendix  presents  a  brief  development  of  the 
mathematical  theory  underlying  the  one-dimensional  GEOFAST 
algorithm.  The  material  is  condensed  from  Ref.  10,  where  fur¬ 
ther  details  may  be  found,  and  is  included  here  for  complete¬ 
ness.  Some  familiarity  with  these  results  is  required  for  an 
understanding  of  the  two-dimensional  theory  discussed  in 
Appendix  B. 


A . 1  TOEPLITZ  AND  CIRCULANT  MATRICES 

Toeplitz  matrices  are  a  class  of  n*n,  real-valued 
matrices  which  includes  the  covariance  matrices  from  station¬ 
ary  random  processes.  More  formally,  a  matrix  T  is  of  the 
Toeplitz  type  if  there  is  a  real-valued  function,  t^ ,  such 
that  the  elements  of  T  obey 

IT] JR  =  tk_j  ,  0  <  j,  k  <  n-1  (A. 1-1) 

A  Toeplitz  matrix  is  defined  by  its  first  row  and  column  since 
[T] =  [T]  ...  If  T  is  an  autocovariance  or  a  cross-co- 
variance  matrix,  then  t£  is  the  corresponding  covariance  func¬ 
tion.  The  matrix  T  is  symmetric  if  the  function  t£  is  even. 

Circulant  matrices  are  a  subset  of  Toeplitz  matrices 
with  an  especially  simple  form  under  Fourier  transformation. 

A  Toeplitz  matrix  is  a  circulant  if  the  function  t£  obeys 
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t-IL  =  tn-£  ’  £  1*  n_1  (A. 1-2) 

A  circulant  matrix  is  defined  by  its  first  row  since  each  row 
is  equal  to  the  row  preceding  it  shifted  one  element  to  the 
right,  with  the  last  element  "wrapped  around"  to  the  first 
place.  For  example,  a  circulant  matrix  with  the  first  row 
(Cq,  c-p  ...,  cn_^)  necessarily  has  its  second  row  given  by 
(cn_^,  Cq,  c^ ,  ...,  cn.2^-  While  circulant  matrices  them¬ 
selves  do  not  appear  frequently,  their  properties  are  useful 
in  dealing  with  the  more  general  Toeplitz  matrices. 


Circulant  matrices  have  the  special  property  that 
they  are  diagonalized  by  the  Discrete  Fourier  Transformation 
(DFT).  The  n-point  DFT  matrix  F  is  the  nxn,  symmetric,  uni- 

JL 

tary,  complex-valued  matrix  whose  elements  are 


F] 


jk 


=  —  exp 

V  n 


2ni  jk 
n 


0  <_  j  ,  k  £  n-1 


(A. 1-3) 


and  its  inverse  is  easily  seen  to  be  F' . ■  The  DFT  of  an  n-vector 
of  complex  numbers,  x,  is  defined  by  the  relations  (Ref.  9) 


x' 


F  x. 


x ' 


(A. 1-4) 


The  matrix  similarity  transformation  corresponding  to  Eq . 
A. 1 -4  is  given  by 


C ' 


F  C  F1 


(A. 1-5) 


where  C  is  an  n*n  matrix  of  complex  numbers  and  C'  is  its  rep¬ 
resentation  in  the  Fourier  transform  domain.  The  usefulness 
of  the  circulant  matrix  definition  depends  largely  on  the  fol- 
lowi ng  resu  1 1  . 

-'•'The  symbol  i  denotes  -J - 1  ,  and  the  superscript  ( t  )  denotes  the 
complex  conjugate  transpose. 
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The  Fourier  transformation,  Eq .  A. 1-5,  of  any  circulant 
matrix  C  is  the  diagonal  matrix 

C'  =  diag(c^ ,  ....  (A. 1-6) 

where  c'  is  the  discrete  power  spectrum 

c '  =  Jn  F^  c  (A. 1-7) 

and  c  corresponds  to  the  first  row  of  C.  Moreover,  any  diagonal 
matrix  is  the  Fourier  transform  of  a  circulant  matrix,  defined 
by  Eqs .  A. 1-6  and  A. 1-7.  If  the  matrix  C  is  symmetric  the  in¬ 
verse  transform  F^  may  be  replaced  by  F  in  Eq.  A. 1-7.  A  proof 
of  this  result  may  be  found  in  Refs.  2  and  4.  This  theorem  is 
also  equivalent  to  the  well  known  relationship  between  discrete 
convolutions  and  discrete  Fourier  transforms  proved  in  Ref. 

20. 


Since  circulant  matrices  are  diagonalized  by  the  DFT 
they  can  be  multiplied  or  inverted  in  order  nlog2n  operations. 
Specifically,  if 


y  =  C  x 


(A. 1-8) 


then 

y'  =  C '  x '  =  diag(c' )  x '  (A. 1-9) 

where  x'  and  y*  are  defined  as  in  Eq .  A. 1-4  and  c'  in  Eq . 

A. 1-7.  It  follows  that 

y '.  =  c  •  x  ,  x =  y'./c\  ( A .  1  - 1 0 ) 

Jj  J  J  J  J  J 

so  that  either  y  or  x  can  be  found  with  two  FFT's,  n  scalar 
multiplications,  and  one  inverse  FFT.  The  matrix  multiplication 
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r 


in  Eq .  A. 1-8  is  also  equivalent  to  the  circular  convolution 


n-1 

y  •  =  /  c,  .  x . 

^  k-j  j 
k=o 


(A. 1-11) 


so  that  fast  multiplication  is  the  same  as  fast  convolution. 


A. 2  THE  TRANSFORM  OF  A  WINDOWED  T0EPL1TZ  MATRIX 

This  section  describes  the  calculation  of  the  Fourier 
transform  of  a  "windowed"  Toeplitz  matrix.  An  algorithm  to 
accomplish  this  for  a  single  superdiagonal  band  in  approxi¬ 
mately  nlogn  computer  operations  is  derived  here  from  the 
result  of  Section  A.l.  A  windowed  Toeplitz  matrix  is  one 
which  has  been  multiplied  by  a  diagonal  matrix,  W,  to  form 

U  =  W  T  W  (A. 2-1) 

The  matrix  W  is  defined  by  a  window  function,  w£ ,  such  that 

W  =  diag( Wq ,  Wj ,  ...,  wn_j)  (A. 2-2) 

The  window  function  is  related  to  the  way  in  which  data  is 
processed,  and  is  described  in  the  next  section.  For  a  gen¬ 
eral  discussion  of  window  functions  in  frequency  domain  data 
processing,  see  Refs.  8  and  19.  Note  that  one  choice  for  W  is 
simply  W  =  I  where  I  is  the  n*n  identity  matrix.  Thus,  the 
class  of  windowed  Toeplitz  matrices  includes  all  Toeplitz 
ma t  rices. 

The  DFT  of  the  windowed  Toeplitz  matrix  in  Eq .  A. 2-1 
i  s 


T  1 


F  W  T  W  F1” 


(A. 2-3) 


The  work  of  deriving  a  useful  element-by-element  expression 
for  this  matrix  is  divided  into  two  parts:  first,  obtaining 
an  expression  for  T'  in  terms  of  transforms  of  diagonal  and 
circulant  matrices  and,  second,  applying  the  circulant  theorem 
of  Section  A.l  to  deduce  an  efficient  computational  formula. 


The  matrix  T  is  extended  to  a  2n*2n  circulant  matrix 
¥  by  the  definitions 


0  <  j,  k  <_  2n-l 


(A. 2-4) 


where 


■  j  0 

(C£-2n 


0  <  £  <_  n-1 
£  =  n 

n+1  <  £  <  2n-l 


(A. 2-5) 


and  c_£  =  ^2n-H  ^or  ^  ^  1  2n-l.  The  matrix  ¥  has  the  parti¬ 

tioned  form 


¥  = 


T 

R 


R 

T. 


(A. 2-6) 


where  R  is  a  Toeplitz  matrix  corresponding  to  the  function 


J 


'£  -n 

(t"-1  +  t-n+l) 


'n-1 
t 


£  +  n 


£  >  0 
£  =  0 
£  <  0 


Also  a  2n*2n  extended  version  of  W  is  defined  by 


(A. 2-7) 
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w  = 


w 

0 


0 

0 


(A. 2-8) 


The  windowed  Toeplitz  matrix  WTW  may  now  be  expressed 
in  terms  of  the  diagonal  matrix  W  and  the  circulant  matrix  T  by 


WTW  0 
0  0 


WTW 


(A. 2-9) 


Introducing  the  n*2n  sampling  matrix  S  defined  by 


1  0  .  .  .  0 
0  0  1.  .  0 


S  = 


0  10 


it  may  be  verified  that 


T' 


~  ~  ~  ~  T 
S  F  W  T  W  F'S1 


(A. 2-10) 


(A. 2-11) 


where  F  is  the  2nx2n  Discrete  Fourier  Transformation  defined 
by  Eq .  A. 1-3  with  n  replaced  by  2n.  Equation  A. 2-11  expresses 
T'  as  the  result  of  sampling  every  other  row  and  column  of  a 
matrix  which  is  the  transform  of  a  product  of  diagonal  window 
matrices  and  a  circulant  matrix. 


The  application  of  the  circulant  properties  stated  in 
Section  A.l  produces  the  desired  expression  ior  T'.  Equation 
A. 2-11  may  be  written  as 

T'  =  S  W'T'W'S1  (A. 2-12) 
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(A. 2-13) 


where  W'  and  f'  are  the  transformed  matrices 

W'  =  F  W  F* 

T'  =  F  T  F*  (A. 2-14) 

Using  Eqs.  A. 1-5  and  A. 1-7  with  F  replaced  by  F 
shows  that  the  first  of  these  matrices  is  a  circulant,  with 
the  elements 


( W’ ]jk  =  Qk_.  ,  0  <  j ,  k  <  2n-l 
where  is  given  by  the  transform 

w 
0 


0/ 


fi 


2n- 1 


=  n  =  — 


/2K 


(A. 2-15) 


(A. 2-16) 


and  w  is  the  n-vector  defined  by  the  window  function 

w  =  (wQ,  wlt  w^)1  (A. 2-17) 

It  also  follows  directly  from  Eqs.  A. 1-5  and  A. 1-7 
that  the  matrix  T'  of  Eq .  A. 2-14  is  diagonal,  of  the  form 

T'  =  diag  ( T q ,  Tlf  T 2n_ i )  (A. 2-18) 

T 

where  t  -  <Tq,  Tj,  x2n_l^  xS  defined  by  the  vector  c  of 

Eq .  A. 2-5  via 


T 


(A. 2-19) 


A  useful  element-by-element  formula  for  T'  is  a  con¬ 
sequence  of  Eq .  A. 2-12.  Substitution  of  the  definitions  from 
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Eqs .  A. 2-10,  A. 2-15,  and  A. 2-18,  with  some  simplification, 
produces  the  expression 


IT'  j 


jk 


2n-l 


E 


Q 


£-2j 


fi£-2k  T £ 


j 

(A. 2-20) 


where  is  the  transform  of  the  window  function  (Eq.  A. 2-16) 
and  is  the  power  spectrum  of  T  (Eq.  A. 2-19).  The  sub¬ 
scripts  in  Eq .  A. 2-20  are  interpreted  modulo  2n.  In  partic¬ 
ular  the  elements  of  the  m11^  superdiagonal  band  of  T'  can  be 
written  (for  0  £  j  <  n-m-1) 


[T '  ] 


j  .  j+m 


2n-l 

E 

£  =  0 


p£-2j (m)  T£ 


(A. 2-21) 


where 

Pk<-0  =  ak  Sk.2n,  (A. 2-22) 

Now  Eq .  A. 2-21  can  be  recognized  as  a  convolution,  or  alterna¬ 
tively  as  the  product  of  a  circulant  matrix  and  a  vector  (Sec¬ 
tion  A.l).  As  a  result  the  band  elements  defined  by  Eq .  A. 2-21 
can  be  computed  in  order  nlog2n  operations. 


A. 3  WINDOWS  AND  BANDING 

The  impact  of  the  window  function  on  the  structure  of 
T'  can  be  seen  from  Eqs.  A. 2-21  and  A. 2-22.  The  elements  of 
the  diagonal  band  are  obtained  by  convolving  the  power 


The  symbol  represents  the  complex  conjugate  of  the  quan¬ 
tity  fi£.  Since  W '  is  hermitian,  Q .  =  fi  .  _  ^ . 
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spectrum  t^  with  the  weighting  function  p^(m),  which  consists 
of  the  transformed  data  window  multiplied  by  its  conjugate 
shifted  2m  elements  to  the  right. 

Two  typical  window  functions  and  their  transforms  are 
shown  in  Fig.  A. 3-1.  These  windows  are  special  cases  of  the 
Kaiser  window  function  which  has  proved  very  useful  in  signal 
processing  applications,  and  is  defined  in  Refs.  15  and  19. 
Transforms  of  data  windows  are  characterized  by  a  mainlobe  and 
a  series  of  sidelobes  of  considerably  reduced  levels.  As  can 
be  seen  from  the  figure,  the  disparity  between  mainlobe  and 
sidelobe  levels  becomes  more  pronounced  as  the  width  of  the 
mainlobe  is  increased.  This  behavior  is  reflected  in  the  nu¬ 
merical  structure  of  the  matrix  T’.  The  elements  of  T'  corre¬ 
sponding  to  the  shaded  area  of  Fig.  A. 3-2  cannot  be  neglected 
because  the  mainlobes  of  Ci  ^  and  2m  overlap  for  these  ele¬ 
ments.  However,  the  elements  of  T'  outside  the  shaded  area  of 
the  figure  are  the  product  of  the  sidelobes  of  together 
with  the  power  spectrum  t£.  When  the  sidelobes  are  made  suf¬ 
ficiently  small,  these  elements  may  be  neglected.  The  key 
idea  in  the  band-diagonal  algorithm  is  the  control  of  out- 
of-band  covariance  elements  by  choice  of  a  data  window  with 
sufficiently  small  sidelobes. 

A  final  transformation  must  be  applied  to  form  the 
frequency  domain  covariance  matrix  into  a  band-diagonal  struc¬ 
ture.  Because  it  is  known  that  the  data  must  always  be  real¬ 
valued,  the  sine-cosine  transform  coefficients  are  sufficient 
to  describe  the  data,  rather  than  the  full,  complex  set  of 
Fourier  coefficients.  The  practical  importance  of  this  last 
transformation  is  the  elimination  of  the  elements  in  the  upper 
right  and  lower  left  corners  of  T'  as  shown  in  Fig.  A. 3-2. 
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Data  Windows  and  Transforms 
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Figure  A. 3-2  Structure  of  Covariance  Matrix 

Under  Fourier  Transformation 

The  transformation  to  sine-cosine  coefficients  corre¬ 
sponds  to  modifying  the  complex  transformation  of  Eq .  A. 2-3  to 
the  real  form 

T'  =  H  F  W  T  W  F^H^  (A. 3-1) 

where  H  is  an  n*n  complex  transformation  from  exponential 
Fourier  coefficients  to  sine-cosine  coefficients.  The  ele¬ 
ments  of  H  are  defined  by 
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The  first  n/2+1  rows  of  H  generate  the  cosine  coefficients  of 
the  data,  while  the  last  n/2-1  rows  generate  the  sine  coefficients. 

The  result  of  using  the  modified  transformation  (Eq. 

A. 3-1)  in  the  definition  of  T'  is  shown  in  Fig.  A. 3-3.  The 
matrix  T'  may  now  be  approximated  by  the  real  band-diagonal 
matrix,  D',  neglecting  the  small  elements  due  to  the  sidelobes 
of  Qg .  In  addition  the  elements  of  D'  may  be  partitioned  into 
blocks  of  dimension  n/2+1  and  n/2-1  so  that  the  dimension  of 
the  largest  system  of  equations  to  be  solved  is  approximately 
halved . 


-■>  )  ROWE 


Figure  A . 3- 3 


Structure  of  Covariance  Matrix  Under 
Sine-Cosine  Transformation 


A. 4  APPROXIMATE  SOLUTION  OF  TOEPLITZ  EQUATIONS 


The  solution  of  the  linear  system 


y  =  I  x 


(A. 4-1) 


can  be  obtained  efficiently  using  the  transformation  of  Eq . 
A. 3-1  when  T  is  a  symmetric  Toeplitz  matrix.  If  A  =  HFW  and 


y'  =  A  y 


(A. 4-2) 


then 

l'  =  T'x 1  ,  T'  =  A  T  A1"  (A. 4-3) 

where  T'  is  symmetric  and  has  the  structure  in  Fig.  A. 3-3. 

This  is  a  generalization  of  the  circulant  case  where  T'  is 
diagonal.  An  efficient  aproximate  solution  to  Eq.  A. 4-3  is 
found  by  replacing  the  exact  matrix  T'  by  its  banded  approxi¬ 
mation  D'.  The  special  structure  of  D'  is  exploited  by  employ¬ 
ing  the  band-diagonal  implementation  of  the  Cholesky  decompo¬ 
sition  algorithm  (Ref.  6,  Section  2.3)  for  the  solution  of 

linear  equations.  The  number  of  numerical  operations  required 

2 

is  of  order  (mB+l)  n  for  large  n,  where  is  the  number  of 
superdiagonal  bands. 

The  magnitude  of  the  band  approximation  error  can  be 
expressed  by  the  quantity 


e 


E(6x  6x) 
E ( xTx ) 


1/2 


(A. 4-4) 


which  is  the  relative  rms  error  in  the  solution  x  to  Eq .  A. 4-1 
due  to  replacing  T'  by  D’  in  Eq .  A. 4-3.  The  error  in  the 
solution  due  to  the  neglected  out-of-band  elements  in  T'  is 
controlled  by  adjusting  the  matrix  bandwidth  m^  through  the 
choice  of  the  data  window  function  defined  in  Section  A. 3. 
Bounds  on  the  computational  error  arising  from  this  source 
were  computed  in  Ref.  10  and  show  that  <  10  is  sufficient 
for  most  applications. 
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There  is,  however,  a  penalty  for  the  use  of  highly 
tapered  data  windows  (large  mQ ) .  It  can  be  shown  that  the 

D 

condition  number  K  of  the  matrix  T'  satisfies  (Ref.  6) 

K(T')  =  A  (T 1 )/A  .  (T')  >  w"2  (A. 4-5) 

max  nun  —  o 

where  w  is  the  smallest  window  coefficient  and  A(T')  denotes 
o 

an  eigenvalue  of  the  matrix  T'.  When  K  is  large  the  matrix  is 
ill-conditioned  and  computing  its  inverse  by  the  Cholesky  meth¬ 
od  may  fail  due  to  the  appearance  of  negative  diagonal  elements. 
To  compensate  for  ill-conditioning  the  equations  A. 4-3  are 
modified  by  adding  a  small  constant  6  to  the  diagonal  so  that 

T’  =  A  T  A*  +  6  1  (A. 4-6) 

The  modified  problem  is  well-conditioned  since  A  .  (T')  >  6 
r  man  — 

and  K  <  A  (T')/6.  By  inverting  the  transformation  (Eq.  A. 4-3) 
the  modified  equations  are  found  to  be 

y  =  (T  +  6  W"2)  x  (A. 4-7) 

where  W  is  the  diagonal  window  matrix. 

When  T  is  a  covariance  matrix  Eq .  A. 4-7  may  be  inter¬ 
preted  as  the  addition  of  uncorrelated  measurement  noise  to 
the  model.  In  the  modified  problem  the  quality  of  the  data  is 
gradually  de-emphasized  as  the  ends  of  the  data  interval  are 
approached,  and  points  for  which  the  signal - to-noise  ratio 

sk  =  tcw2/6  (A. 4-8) 

is  less  than  one,  may  be  regarded  as  having  been  lost.  The 
percentage  of  data  points  for  which  is  less  than  one  de¬ 
fines  a  performance  index,  called  "data  de-emphasis,"  which  is 


58 


a  function  of  6  and  mg. 
that  a  relative  error  of 
data  de-emphasis  of  less 


The  results  presented  in  Ref.  10  show 
less  than  1%  may  be  achieved  with  a 
than  10%  for  typical  applications. 


f 

4 


j 
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APPENDIX  B 

GEOFAST:  TWO-DIMENSIONAL  THEORY 


This  appendix  extends  the  mathematical  theory  under¬ 
lying  the  GEOFAST  algorithm  to  two  dimensions.  Since  a  basic 
understanding  of  the  one-dimensional  theory  (Ref.  10)  is  re¬ 
quired,  the  necessary  results  are  developed  in  Appendix  A. 
Thus  the  two  appendices  taken  together  are  self-contained. 


B.l  BLOCK  CIRCULANT  MATRICES  AND  KRONECKER  PRODUCTS 

The  relationship  between  circulant  matrices  and  the 

Discrete  Fourier  Transform  (DFT)  is  described  in  Appendix  A. 

“V 

The  DFT  matrix  in  one  dimension  is  defined  by 


[FJ,,  =  wjk 


w  =  exp 


(B.l-1) 


where  u»  is  the  complex  nL  root  of  unity  and  n  is  the  number 
of  points.  A  circulant  matrix  satisfies 


C  ^  j  k  c k- j (mod  n ) 


0  £  j,k  <  n-1 


(B.l-2) 


and  is  defined  by  its  first  row  "circulated".  The  significance 
of  these  definitions  is  that  any  circulant  matrix  is  diagonal¬ 
ized  by  the  DFT,  that  is. 


C'  =  F  C  F^  =  /n  diag  |F^c] 


(B.l-3) 


■The  symbol  i  denotes  /- 1  and  the  superscript  (t)  denotes  the 
complex  conjugate  transpose. 
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is  a  diagonal  matrix  (Refs.  2  and  4).  This  relationship  ex¬ 
tends  to  two  dimensions  with  the  appropriate  definitions  as 
shown  below. 


B.1.1  The  DFT  as  a  Kronecker  Product 


The  DFT  in  two  dimensions  is  defined  by 


nl~l  n2_l 


pj  qk 

pq  ^  Jk  “1  “2 

vnln2  j=o  k=o 


(B.l-4) 


where  Zj^  is  an  n^*n2  data  matrix  and  denote  the  com¬ 

plex  roots  of  unity 


2n  i 

2n  i 

=  exp 

nl_ 

,  u)2  =  exp 

n2 

(B.l-5) 


Analogous  to  the  DFT  matrix  in  Eq.  B.l-1,  the  operator  defined 

by  Eq.  B.l-4  is  a  tensor  of  fourth  rank  whose  components  are 

of  the  form  a^  . 

pqjk 

It  is  more  convenient  to  represent  the  DFT  in  matrix 
form  by  expressing  z^  (row-wise)  as  a  vector  z,  such  that 


=  F  z 


(B.l-6) 


is  equivalent  to  Eq.  B.l-4.  When  this  is  done  the  n^n2*n^n2 
matrix  F  is  seen  to  assume  a  simple  block  structure,  namely. 


fl ,11  F2  ! 

f 1 , 12  F2 

f 

1 

1 

1 

f 1 ,21  F2  ! 

„ _ J. 

f 1 , 22  F2 

1 

1 

1 

L _ 

(B.l-7) 

1 

1 

1 

1 

1 

1 
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where  F2  is  the  ^-point  DFT  matrix  defined  by 


[F,l 


=  f. 


2  J  j  k  ~  2  ,  j  k 


1 


UJ 


jk 


(B.  1-8) 


and  similarly  for  F^.  Each  block  in  Eq .  B.l-7  is  of  size  n2*ri2 
and  there  are  n^  blocks  in  each  row  and  column.  This  special 
structure  is  denoted  by 


F  =  F]_  ®  F2  (B.l-9) 

and  termed  the  Kronecker  product  of  F^  and  F2  (Ref.  2).  Kronecker 
products  define  a  special  class  of  tensors  with  the  important 
properties 


(A®  B)  (C  ®  D)  =  (AC  ®  BD) 

(A®  B)1  =  (A1  ®  Bf )  (B.l-10) 

(A®  B)'1  =  (A-1  ®  B"1) 

where  the  matrices  A,  C  and  B,  D  are  assumed  to  be  conformably 
dimensioned . 

For  computational  purposes  it  is  convenient  to  eval¬ 
uate  the  multiplication  of  a  vector  by  a  Kronecker  product 
matrix  in  terms  of  ordinary  matrix  operations.  The  two  prod¬ 
uct  forms 


y  =  (A  ®  B )  x 
Y  =  A  X  BT 


(B.l-11) 


are  equivalent  if  x,  y  are  the  row-by-row  vector  representations 
of  the  matrices  X,  Y.  It  is  assumed  that  A  is  n^xn^  ,  B  is 
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i 

« 


n2Xn2’  an<^  ^  are  n]^Xn2*  This  equivalence  can  be  derived 
from  the  corresponding  element-by-element  formulas  such  as  Eq. 
B.l-4. 


B.1.2  Diagonalization  of  Block  Circulant  Matrices 


By  a  block  circulant  matrix  (with  circulant  blocks) 
is  meant  an  n^^xn-j^  matrix  of  the  form 


C  = 


-2 


♦ 


(B. 1-12) 


where  each  of  the  n^  block  matrices  is  a  circulant  of  dimen 
sion  n2*n2'  A  special  case  of  a  block  circulant  matrix  is  the 
Kronecker  product  of  two  circulants  ,  C2  of  dimension  n^ ,  n2 


C  =  cx  ®  c2 


(B.l-13) 


In  this  case  Eqs.  B.l-3,  B.l-9  and  B.l-10  show  that 

C'  =  F  C  Ft  =  (F]_  C:  F^)®  (F2  C2  F2t)  (B.l-14) 

is  a  diagonal  matrix.  The  diagonal  entries  of  C'  are  of  the 
products  of  the  diagonal  entries  of  Cj  ,  . 


Of  more  interest  is  the  fact  that  C'  is  diagonal  for 
any  block  circulant  C.  This  theorem  is  derived  below.  Let 
denote  the  primitive  circulant  matrix  of  dimension 


1 


4 


1 

j 
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P1  -- 


0  0 
0  0 


(B. 1-15) 


T 

The  matrix  is  orthogonal  (P^P^  =  ^  an<^  generates  a  cyclic 
group  of  order  n-^ 


!Pj,  ?l 


n-,  - 1  n,  n 

P1  1  ,  P1  1  =  P^  =  1}  (B . 1-16) 


since  the  kL  power  of  P,  is  a  circulant  whose  first  row  is 
til 

the  kL  unit  vector.  Furthermore  any  block  circulant  matrix 
may  be  represented  as  the  sum 


"r1 


C  =  £  (Pik©ck 


(B.l-17) 


where  each  is  an  n2xn2  circulant  (see  Eq .  B.l-12).  The 
theorem  follows  since 


"r1 

=  F  C  FT  =  (Fj  Pjk  F^)©  (F2  Ck  F2t)  (B.l-18) 


and  each  factor  in  parentheses  is  the  DFT  of  a  circulant 


An  explicit  expression  for  the  diagonal  entries  of  C' 

K. 

can  be  obtained  by  reference  to  Eq  .  B .  1  -  'i  .  Since  P,  is  a 

t  h 

circulant  whose  first  row  is  the  kL  unit  vector,  e^,  and  C^, 
is  a  circulant  with  first  row  c^,  it  follows  that 
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(B. 1-19) 


F1  P1  F1  =  diaglvrn1  F^  — k  ^ 

F i  ^k  F ^  -  diag  (  ^2  ^2  — -k^ 

Eq.  B.l-18  then  becomes 

C'  =  diag  (  (-/n^  F^  ®  </n2  F2)^  c] 

where  c  is  the  first  row  of  C,  or  equivalently 

nrl 

c  =  2  -k  ®  -k 

k-o 


(B.l-20) 


(B.l-21) 


Note  that  these  two-dimensional  relations  reduce  to  the  one¬ 
dimensional  forms  if  n2  =  1,  and  also  may  be  readily  extended 
to  higher  dimensions  in  a  natural  way. 

The  converse  theorem  is  also  needed,  namely  that  the 
DFT  of  a  diagonal  matrix  D  is  a  block  circulant.  Since  any 
nln2~SC*uare  diagonal  matrix  can  be  represented  as 

"r1 

D=^  [diag(€k)l  ®  [diag  (dk)l  (B.l-22) 

k=o 


where  e^  is  a  unit  vector  of  dimension  n^  and  is  a  vector 
of  dimension  n9 , 

nl_1 

C  =  Ff  D  F  =  ^  [y\  diag(ek)  Fj ]  ®  (F*  diag(dk>  F2 ] 

k=°  (B.l-23) 

The  matrices  in  square  brackets  are  circulants  by  virtue  of 
the  one-dimensional  theory  (Appendix  A).  But  the  Kronecker 
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product  of  two  circulants  is  block  circulant,  and  this  prop¬ 
erty  is  preserved  under  addition.  Since  C  is  a  block  circu¬ 
lant  Eq .  B.l-20  can  be  inverted  to  obtain 

c  =(  —  Fi  ®  —  )  d  (B.l-24) 

where  D  =  diag(d) .  In  particular  if  d  :  0  d2  then 

C  =  C,  ®  C2  with  the  obvious  definitions. 


Since  block  circulants  are  diagonalized  by  the  DFT 
they  can  be  multiplied  or  inverted  in  order  Nlog2N  operations 
in  complete  analogy  with  the  one-dimensional  case  (Appendix  A) . 
If 


then 


y  =  C  x 


(B.l-25) 


y '  =  diag( c ' )  x ' 


(B . 1-26) 


where  prime 
Thus  either 
plications , 


denotes  transform  and  c'  is  given  by  Eq.  B.l-20. 
x  or  y  can  be  found  with  two  FFT's,  N  scalar  multi- 
and  one  inverse  FFT. 


The  matrix  multiplication  in  Eq .  B.l-25  is  also  equiv¬ 
alent  to  a  two-dimensional  circular  convolution 


jk 


n j -  1  n2 -  1 

£  £ 

p=o  q=o 


c  ,  x 

p-j.q-k  pq 


(B. 1-27 ) 


where  X,  Y  are  the  matrix  forms  of  x,  y, 
form  the  vectors 


and  the  elements  c.. 

Jk 
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* 


(B.l-28) 


which  define  the  circulants  in  Eqs.  B.l-12.  In  the  special 
case  that  C  =  A®  B,  with  A  and  B  circulant,  Eq .  B.l-27  is 
simplified  to 

n^-1  n2"l 

yjk  =  2  Vj  £  bq-k  xpq  (B.  1-29) 

p=o  q=o 

which  is  a  sequence  of  row  and  column  convolutions  (see  Eq . 
B.l-11) . 


B . 2  THE  TRANSFORM  OF  A  WINDOWED  BLOCK  TOEPLITZ  MATRIX 

B.2.1  Matrix  Representation  of  the  Transform 

A  block  Toeplitz  matrix  (with  Toeplitz  blocks)  is  an 
n-j^xn^j^  matrix  of  the  form  (see  Eq.  B.l-12) 


1 

— 1 

T0  i  T1 

T»r> 

-  - 

T-1  i  T0 

V* 

i 

-  - 

1 

_  _ 

T-(nrl)  !  T-(nr2) 

T0 

—  1 

where  each  block  T^  is  a  Toeplitz  matrix  of  dimension  n2xri2' 
Such  matrices  arise  naturally  as  the  covariance  matrix  of  a 
stationary  process  on  a  planar  rectangular  grid.  A  special 
case,  corresponding  to  a  separable  covariance  function,  is  the 
Kronecker  product 
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T  =  T1  0  T2 


( B . 2-2  ) 


where  ,  T2  are  Toeplitz  of  dmension  ,  n 2  * 

Toeplitz  matrices  in  one  dimension  are  discussed  in 
Appendix  A,  in  particular  the  transform  of  a  windowed  n2xn2 
Toepl i tz  matrix 


T 2 '  =  F2  W2  T2  W2  F2' 


(B.2-3) 


where  \v'2  is  a  diagonal  window  matrix. 

The  development  given  in  Section  A. 2  extends  directly 
to  two  dimensions.  Define  the  n^n2*n^n2  diagonal  window  matrix 
W  as  the  Kronecker  product  of  two  one-dimensional  windows 


W  =  Wx  ®  w2 


(B.2-4) 


This  is  a  natural  choice  of  window  function  on  a  rectangular 
grid.  The  transform  of  a  windowed  block  Toeplitz  matrix  anal¬ 
ogous  to  Eq .  B.2-3  is  the  n^n2xn^n2  matrix 


T'  =  F  W  T  W  F' 


(B.2-5) 


where  T  is  defined  in  Eq .  B.2-1,  W  in  Eq .  B.2-4  and  F  in  Eq . 

B.  1-9. 

Extend  each  Toeplitz  block  T^  in  Eq .  B.2-1  to  a  2n2x2n2 
circulant  in  the  manner  of  Appendix  A 


A  rk1 


L  Rk  TrJ 


(B.2-6) 
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If  Che  Toeplitz  matrix  is  defined  by  the  vector 

t=  ...  ,  t_1,  t0,  tj,  ...  ,  tfl2_ i_  1  (B.2-7) 


then  R^  is  the  complementary  Toeplitz  matrix  defined  by 


-  '  lt-l*  •••’  t-(n2-l)’  I(t-(n2-l)  +  tn2-l)  ’  ^-l’  C1 ^ 

(B .2-8) 


Now  extend  T  to  a  An^n2x4n^n?  block  circulant  in  the  same  manner 


T  = 


T  R 

a-  -J 

R  '  T 


(B.2-9) 


where  T  is  the  2njn2*2njn2  block  Toeplitz  matrix 

1 


0 


T 


1 


nr] 


T 


T 


-1 


(B.2-10) 


-(n, -1  ) 


L  _(ni 


and  R  is  the  complementary  block  Toeplitz  matrix 

1 


R 


i(T  .  .  ♦T  ,  )  T  . 

2  - ( n  j  -  1  )  n  j -  1  n ,  -  1 


T1  ! 


T 


( n  i  -  1  ) 


J 


(B.2-11) 
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The  extension  of  the  window  matrix  and  DFT  to  compat¬ 
ible  size  (4n^n2)  is  straightforward.  That  is, 

W  0' 

W  =  =  W,  ®  W9  (B.2-12) 

0  0  J  1  z 

where  ,  W9  are  one-dimensional  extensions  to  size  2n^  ,  2n2 , 
and 

F  =  F:  ©  F2  (B.2-13) 

where  F^ ,  are  DFT ' s  of  size  2n^  ,  2n2 .  Thus  Eq .  B.2-5  can 
be  written 

T'  =SFWTWFST  (B.2-14) 

where  S  is  a  two-dimensional  sampling  matrix 

S  =  S1  ®  S2  (B.2-15) 

formed  from  two  sampling  matrices  as  defined  in  Appendix  A. 
Finally, 

T'  =  S  W'T'W'S7  (B. 2-16) 

where  the  transformed  matrices 

W*  =  F  W  Ft  ,  T’  =  F  T  F1  (B.2-17) 


are  block  circulant  and  diagonal,  respectively,  by  the  theorems 
o f  Sec  t ion  B . 1 . 


B.2.2  Efficient  Calculation  of  the  Transform 


Equation  B.2-16  is  the  direct  analog  of  the  one-dimen¬ 
sional  Eq.  A. 2-12,  and  can  be  used  to  exhibit  an  element-by¬ 
element  formula  analogous  to  Eq .  A. 2-20.  It  is  instructive  to 
express  Eq.  B.2-16  in  terms  of  the  underlying  one-dimensional 
transforms  so  that  the  one-dimensional  theory  can  be  applied. 
This  is  especially  easy  when  the  covariance  function  is  separ¬ 
able,  since  then  Eqs .  B.2-2  and  B.2-12  to  B.2-15  show  that  the 
computation  "factors"  completely  to 

T’  =  T1'  ®  T2'  (B.2-18) 

where 

Tl'  =  S1  VW8!1  (B .  2-19) 

is  the  one-dimensional  transform  of  T^ ,  and  similarly  for  T^. 
Thus  only  two  one-dimensional  transforms  are  required  in  this 
case,  and  their  banded  approximations  may  be  calculated  as 
described  in  Appendix  A  in  the  order  of  mg(n^log2n-^  +  n2log2n2) 
operations . 

In  the  non-separable  case  this  factorization  is  only 
partial  in  the  sense  that  T'  in  Eq .  B.2-16  does  not  factor  into 
a  Kronecker  product.  However,  T'  is  a  diagonal  matrix  and  can 
be  represented  as  a  sum  of  products  of  diagonal  matrices 

2n^-l  2n2~l 

T'  =  [diag(ep)  ®  diag(e^)  )  Tpfj  (B.2-20) 

p=o  q=o 

x 

This  assumes  that  both  windows  ,  W9  have  the  same  bandwidth, 
mn,  which  is  not  required. 
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f  Vi  t  Vi 

where  e  ,  e^  are  p  ,  q  unit  vectors  of  dimension  2n^ ,  2n2 , 
and  t  are  the  elements  of  a  2n,x2n0  matrix  containing  the 
diagonal  entries  of  T'.  Specifically,  let  t_  be  the  4nin2~ 
vector  corresponding  row-wise  to  Then  from  Eq .  B.l-20 


t  = 


Fx  ®  V2K2  f2 


(B .2-21) 


where  c  is  the  first  row  of  the  circulant  T  constructed  in 
Eqs.  B.2-6  to  B.2-11.  The  vector  c  can  be  defined  explicitly 
in  terms  of  the  first  row  and  column  of  T  by  Eq.  B.2-1. 

With  the  representation  of  Eq .  B.2-20  the  transform 
Eq .  B.2-16  becomes 


2n^-l  2n2~l 


=  E  E  <Ap®V 


(B.2-22) 


p=o  q=o 


where 


A’  =  S1  W]  diag(e  )  W] 


B-  =  S2  W’  diag( e  )  sj 


(B .2-23) 


The  calculation  of  A'  and  B'  is  quite  simple  since  diag(e  )  =  e 

P  4  F  6  -p  -p 

imp  1 i  es  that 


A'  = 
P 


a  a1 
~P  “P 


a  =  S,  W’  e 
-p  1  1  -p 


(B.2-24) 


Note  that  Wj  -  Fj  is  herinitian  so  that  (Wj  =  Wj 
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Furthermore,  the  matrix  Wj  is  a  circulant  whose  first  row  is 
given  by 


n,  =  — —  fJ 
_1  1 


*1 


(B.2-25) 


so  that  the  vectors  a^  can  all  be  computed  with  a  single  FFT. 
Explicitly , 


-jp  =  '*il2j,p  =  nl,p-2j  <B-2-26> 

where  negative  subscripts  are  interpreted  modulo  2n-^ ,  and 


"k  ajp  akp  fil,p-2j  nl,p-2k 


(B. 2-27 ) 


where  ft,  .  is  the  complex  conjugate  of  ft,  . . 

i  .  J  1 » J 

Since  ft-^  is  the  transform  of  the  window  w^  the  matrix 

Ap  is  approximately  banded  with  bandwidth  trig  (see  Appendix  A). 

The  representation  of  Eq .  B.2-22  then  shows  that  the  matrix  T' 

is  approximated  by  a  block-banded  matrix  (with  banded  blocks) 

when  Ap  and  are  replaced  by  their  (one-dimensional)  banded 

approximations.  This  follows  from  the  fact  that  the  Kronecker 

product  of  banded  matrices  is  block-banded  (see  Eq .  B.l-7), 

and  that  this  property  is  preserved  under  addition.  The  total 

number  of  non-zero  elements  in  the  block-banded  form  of  T'  is 

2 

less  than  (2irig+l)  N,  and  for  each  row  j  (0  <  j  £  N-l)  the  band 
elements  occur  in  the  columns  j±m  where  m  =  m-j^  +  , 

n»2  nig ) . 

The  elements  of  T'  may  be  efficiently  computed  bandwis 
as  in  the  one-dimensional  case.  In  fact  if  j  =  j^r^  +  j'2  t*len 


2n^-l  2n2~l 

IT’)-  =  V  V'  a  (m,  )  (3  0.  (m0)  x 

J  »J+m  Z-rf  Z-,  P"2j ,  1'  ^q-2j2v  2  pq 

p=o  q=o 

(B.2-28) 

where  0  <  <  n-^-1  »  0  1  J2  £  n2”^’  anc* 

ak(m)  =  °l,k  Ql,k-2m  (B  .  2-29  > 

pk(m)  =  Q2>k  ft2,k-2m 

The  analogy  with  the  one-dimensional  formula  is  completed  by 
defining 


p(m)  =  a (m1 )  ©  £  (m2> 


(B. 2-30) 


so  that 


[T'lj 


,  j+m 


4nin2-l 

2  P£-2j(m)  T£ 

£=o 


(B.2-31) 


where  T_  is  defined  by  Eq .  (B.2-21).  Equation  B.2-28  can  be 

recognized  as  a  convolution  or  circulant  multiplication  (see 

Eq .  B.l-29)  and  can  therefore  be  computed  in  the  order  of  Nlog2N 

operations.  The  total  workload  to  compute  all  bands  in  the 

2 

approximation  is  therefore  of  order  m^  Nlog2N. 


In  practice  the  complex  transforms  implied  by  Eqs. 
B.2-25  and  B.2-29  are  replaced  by  real  transforms.  This  amounts 
to  replacing  the  complex  DFT  matrix  F  in  Eq.  B.2-5  by  the  real 
OFT  matrix  HF  where 


H  =  Hl  ®  H2 


(B.2-32) 
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is  a  complex  n^r^xn^r^  matrix,  and  H-p  H2  are  equal  to  the 
matrix  H  defined  in  Appendix  A  (Eq.  A. 3-2)  for  n  =  n,  ,  n0 . 

+  +  +  L  L. 

Similarly,  F  is  replaced  by  F 7 H '  so  that  Eq.  B.2-5  is  formal¬ 
ly  identical  to  Eq.  A. 3-1.  The  result  is  that  the  banded  blocks 
of  T'  have  the  structure  described  for  the  one-dimensional 
case  in  Appendix  A. 


B . 3  APPROXIMATE  SOLUTION  OF  BLOCK  TOEPLITZ  EQUATIONS 
B.3.1  Fast  Multiplication 

There  are  two  possible  ways  to  obtain  a  "fast"  multi¬ 
plication  algorithm  for  block  Toeplitz  matrices,  both  of  which 
are  of  order  Nlog2N  in  complexity.  The  simplest  method  is  the 
application  of  fast  convolution  (Section  B.1.2)  to  the  circulant 
matrix  T  (Eq.  B.2-9)  as  described  in  Ref.  10  for  one  dimension. 

This  method  is  rejected  here  since  experience  with  the  one¬ 
dimensional  algorithm  has  shown  that  estimation  errors  are 
poorly  behaved  with  this  technique.  The  reason  is  that  the 

consistency  in  the  statistical  modelling  (Ref.  13),  which  im-  i 

plicitly  defines  the  relation  between  the  estimates  and  the 
data,  is  destroyed  if  the  auto-covariance  is  band-approximated 
(see  Section  B.3.2)  and  the  cross-covariance  is  not. 

- 

Therefore,  it  is  assumed  that  the  matrix  to  be  multi¬ 
plied  is  replaced  by  its  block-banded  approximation.  A  straight¬ 
forward  block-banded  multiplication  algorithm  then  requires  on  ; 

2  I 

the  order  of  ntgN  operations,  in  addition  to  computing  the  approx¬ 
imation,  since  the  total  number  of  non-zero  elements  in  a  row 

2  2 
is  nig.  The  workload  is  thus  still  dominated  by  nig  Nlog2N. 
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B.3.2  Fast  Inversion  of  Separable  Matrices 

Consider  the  linear  system 

2  =  T  x  (B.3-1) 

where  T  is  block  Toeplitz  of  dimension  0^02X0^2.  Its  fre¬ 
quency  domain  counterpart  is 

2*  =  T'x '  (B.3-2) 

where 

r 

T' 


=  A  2 
=  A  T 


x  =  A^x ' 


(B. 3-3) 


and  A  is  the  transformation  matrix  defined  by 

A  =  H  F  W  (B.3-4) 


where  H,  F,  and  W  are  defined  in  Section  B.2.  The  solution  of 
Eq .  B.3-1  for  x  reduces  to  the  solution  of  Eq .  B.3-2  for  x' . 


If  T  is  separable  and  symmetric  then  from  Eq .  B.2-18 


2'  =  (Tj  ©  T£)  x' 


(B.3-5) 


or  using  Eqs.  B.l-10  and  B.l-11 


x '  =  [ (Tj )'J  ®  (T^)"1 ]  2' 


X’  =  (T j  ) ” 1  V ' ( T 2  )  ~ 1 


(B.3-6) 

(B.3-7) 


where  x',  2'  are  the  row-wise  elements  of  the  nj*n2  matrices 
X'.  V.  Equation  B.3-7  can  be  further  expanded  to  two  simple 
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matrix  inversions 


V'  =  (Ti )  1 ( Y ' )T 


X'  =  (Tj)"1(V')T 


(B.3-8) 


Since  Tj  ,  are  one-dimensional  transforms  of  Toeplitz 
matrices  they  may  be  replaced  by  their  banded  approximations 
and  inverted  by  the  Cholesky  method  as  described  in  Appendix 
A.  However,  the  right-hand-sides  of  Eq .  B.3-8  are  now  matrices, 
resulting  in  a  workload  for  both  equations  of  order 


q  =  mBn2(mB+n1)  +  mBn1(mB+n2)  (B.3-9) 

where  mB  is  the  bandwidth  of  the  approximation.  If  n^  =  n2 
then  q  is  of  order  2mgN. 

B.3.3  Iterative  Solution  for  Nonseparable  Matrices 


When  the  matrix  T  is  not  separable  the  transform  T' 
cannot  be  expressed  as  a  Kronecker  product.  However,  there 
exists  a  natural  approximation  to  T  which  is  separable,  and 
may  be  used  as  the  basis  for  an  iterative  technique.  If  T  is 
a  symmetric  block  Toeplitz  matrix  (Eq.  B.2-1)  it  is  defined  by 
the  sequence  of  vectors 


-k  _  [t-(n2-l),k 


'ok 


(B.3-10) 


which  in  turn  define  the  blocks  for  k  =  0,  1,  ...,  n^-1. 

The  matrix  T  is  separable  if 

ik  =  sk£o  (B.3-11) 
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since  then  T  =  SQ  ®  Tq  where  SQ  is  a  symmetric  Toeplitz  matrix 
defined  by  the  vector 

sQ  =  11, Sj,  ...  ,  sni_iJ  (B.3-12) 

A  natural  Kronecker  approximation  is  obtained  by  choosing  s^  = 
tok//too‘  With  this  choice  the  approximation  agrees  with  T 
identically  for  the  main  diagonal  block  T  ,  and  for  the  main 
diagonal  within  each  block  T^.  The  error  in  the  approximation 
is  determined  by  Eq.  B.3-11. 

Let  D  =  so  (x)  Tq  and  express  T  as 

T  =  D  +  E 

where  E  is  the  error  in  the  separable  matrix  approximation. 
Since  T  and  D  are  Toeplitz  so  is  E,  being  defined  by  the 
vectors 


Sk  =  ik  ■  skio 


(B.3-13) 


The  transformed  linear  system  (Eq.  B.3-2)  becomes 


y'  =  (D’  +  E'  )  x'  (B. 3-14) 

where  O'  is  the  Kronecker  product  ®  T^ ,  and  E'  =  T'  -  D'. 

An  Iterative  solution  to  Eq.  B.3-14  can  be  based  on 
the  identity 

x'  =  (D ' )'1 ly'  -  E ' x '  ]  (B.3-15) 

which  leads  to  the  recursion 
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x'(k+l)  =  ( D ’ )_1 ly '  -  E’  x ' (k) ]  (B. 3-16) 

where  k  is  the  iteration  step  number  and  x'(0)  is  any  conven¬ 
ient  initial  solution  (perhaps  the  zero  vector).  Since  D'  may 

be  inverted  in  order  mnN  operations  from  Section  B.3.2,  and 

"  2 

the  multiplication  by  E'  carried  out  in  order  mRN  from  Section 

d  2 

B.3.1,  the  workload  represented  by  Eq .  B.3-16  is  of  order  msmgN 
where  mg  is  the  number  of  iteration  steps. 

Iterative  methods  for  linear  systems  are  discussed  in 
Ref.  12.  The  recursion,  Eq .  B.3-16,  converges  to  the  solution 
x'  of  Eq.  B.3-14  if  the  approximation  D*  is  sufficiently  close 
to  T’ .  Let 


6x ' (k)  =  x ' (k)-x 1 


(B.3-17) 


4-  y. 

be  the  error  in  the  k  iterate.  Then 


fix’ (k)  =  Mk  fix' (0)  ,  M  =  (D ' ) _1E ' 


(B.3-18) 


and  the  iteration  converges  if  any  matrix  norm  of  M  is  less 

than  one,  in  particular  if  the  eigenvalues  of  M  are  all  less 

than  one  in  absolute  value.  An  expression  for  the  convergence 

rate  in  terms  of  the  number  of  steps,  m  ,  required  to  reduce 

-  B  ^ 

the  initial  error  by  10  p  is  given  by 


m 


s 


P 

-losiop 


,  p  =  max 
i 


(B.3-19) 


where  A^{M)  denotes  the  eigenvalues  of  the  matrix  M.  For  a 
useful  algorithm  mg  should  be  less  than  10  when  p  =  2. 


Since  both  T'  and  D’  are  normally  symmetric  and  posi¬ 
tive  definite,  a  simpler  convergence  criterion  is  (Ref.  12)  that 
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(2D' -T' )  =  D' ( I-M) 


(B.3-20) 


is  a  positive  definite  matrix.  This  holds  when  the  Toeplitz 
matrix  (2D-T)  is  itself  positive  definite,  or  in  particular  if 

max|A^(T)|<  2  min|A.(S  )[•  min|A.(T  )l  (B.3-21) 

i  i  *  1  i  1  1 

since  D  is  a  Kronecker  product.  In  this  case  the  rate  of  con¬ 
vergence  can  be  estimated  using  Eq .  B.3-19  and  the  inequalities 

min  |  Aj(T)|  /max  |  A  ^  ( D )  j  <  1  +  A  ^  (M)  <  max  j  A  ^  (T)  j  /min  J  A  ^  ( D ) 

(B.3-22) 

The  convergence  of  the  iteration  can  be  guaranteed, 

and  an  optimum  rate  of  convergence  obtained,  by  introducing  a 

relaxation  parameter,  a.  Let  D  =  a  where  D  =  ©  T  and 

o  o  o  o 

a  is  a  positive  real  number.  Then  (2D-T)  is  a  positive-definite 
matrix  for  sufficiently  small  a  so  that  convergence  is  assured. 
Substitution  shows  that  Eq .  B.3-16  can  be  rewritten 

x  '  (  k+1 )  =  x  '  (  k )  +  a(DM_1  (y  -  T'x'(k)]  (B.3-23) 

which  is  in  the  form  of  a  relaxation  iteration. 

Analysis  of  the  convergence  rate  of  Eq .  B.3-23  shows 
that  the  optimum  choice  of  a  is  given  by 

«c>  =  2(qj  +  q2 ) (B.3-2A) 

where  q,  and  q,;  are  the  smallest  and  largest  eigenvalues  of 

_  i  1 

<D  T).  The  q's  may  be  estimated  from  upper  and  lower  bounds 

T  T 

for  the  real,  positive  ratio  xlTx/x1D  x.  The  convergence  fac¬ 
tor  corresponding  to  Eq .  B.3-24  is  given  by 
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po  =  ( ^2  "  +  (B  .  3-25 ) 

which  is  less  than  one  for  any  T,  and  reaches  zero  in  the 
separable  case  T  =  Dq.  The  number  of  steps  required  to  reduce 
the  error  by  a  factor  of  100  (see  Eq .  B.3-19)  is  less  than  ten 
when  the  condition  number  of  (D  ^T),  given  by  is  less 

than  about  five.  Furthermore,  it  can  be  shown  that  the  number 
of  steps  increases  linearly  with  the  condition  number. 
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